#pragma once

// Global variables
// Note the maxim number of elements of the array 
double ZETA_GLB[1000];
double ZETACLUST_GLB[1000][1000];
double X_GLB[1000];
int NCLUST_GLB;
int NX_GLB;

double CHI2_lin_chiralfit_GP(double A, double B, double C) {
        double CHI2=0;
        double F1;
        int NPAR=0;
	double *clust_aux=malloc(sizeof(double)*(NCLUST_GLB+1));
	for (int im=0; im<NX_GLB; im++)
	{
                NPAR=NPAR+1;
                for (int iclust=0; iclust<NCLUST_GLB+1; iclust++) clust_aux[iclust]=ZETACLUST_GLB[im][iclust];
                F1=(ZETA_GLB[im]-(A+B*X_GLB[im]+C/X_GLB[im]))/error_jack_real(clust_aux,NCLUST_GLB);
                CHI2=CHI2+pow(F1,2);
        }
	return CHI2/(NPAR-3);
//return CHI2;
}

double CHI2_lin_chiralfit_2GP(double A, double B, double C, double D) {
        double CHI2=0;
        double F1;
        int NPAR=0;
        double *clust_aux=malloc(sizeof(double)*(NCLUST_GLB+1));
        for (int im=0; im<NX_GLB; im++)
        {
                NPAR=NPAR+1;
                for (int iclust=0; iclust<NCLUST_GLB+1; iclust++) clust_aux[iclust]=ZETACLUST_GLB[im][iclust];
                F1=(ZETA_GLB[im]-(A+B*X_GLB[im]+C/X_GLB[im]+D/pow(X_GLB[im],2)))/error_jack_real(clust_aux,NCLUST_GLB);
                CHI2=CHI2+pow(F1,2);
        }
        if ((NPAR-4)==0) NPAR=NPAR+1; //Avoid divide by 0
        return CHI2/(NPAR-4);
//return CHI2;
}

void minuitFunction_lin_chiralfit_GP(int &nDim, double *gout, double &result, double *par, int flg) {
result = CHI2_lin_chiralfit_GP(par[0], par[1] , par[2]);
}
void minuitFunction_lin_chiralfit_2GP(int &nDim, double *gout, double &result, double *par, int flg) {
result = CHI2_lin_chiralfit_2GP(par[0], par[1] , par[2], par[3]);
}



void lin_chiralfit_withGP(double *Z, double **Z_CLUST, double *X, int NX, int Nclust,double &BestA,double &BestB,double &BestC,double &chi2)
{

//Fit to ZETA=A+B*MPS^2+C/MPS^2
//MPS^2 propto (m1+m2)
// M=m1+m2 if POL
// M=MPS if MPS

// Reorganize data in global variables
NX_GLB=NX;
NCLUST_GLB=Nclust;
for (int ix=0; ix<NX_GLB; ix++)
{
	ZETA_GLB[ix]=Z[ix];
	X_GLB[ix]=X[ix];
	for (int iclust=0; iclust<NCLUST_GLB+1; iclust++) ZETACLUST_GLB[ix][iclust]=Z_CLUST[ix][iclust];
}

// Minimize

TFitter* minimizer = new TFitter(2);
// MAKE IT QUIET!!

double p1 = -1;
minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);

// Tell the minimizer about the function to be minimzed
minimizer->SetFCN(minuitFunction_lin_chiralfit_GP);

// Define the parameters
//arg1 – parameter number
//arg2 – parameter name
//arg3 – first guess at parameter value
//arg4 – estimated distance to minimum
//arg5, arg6 – ignore for now
minimizer->SetParameter(0,"A",1.5,0.001,0,0);
minimizer->SetParameter(1,"B",0.0,0.001,0.0,0);
minimizer->SetParameter(2,"C",0.0,0.001,0.0,0);


// Run the simplex minimizer to get close to the minimum
minimizer->ExecuteCommand("SIMPLEX",0,0);
// Run the migrad minimizer (an extended Powell's method) to improve the
// fit.
minimizer->ExecuteCommand("MIGRAD",0,0);
// Get the best fit values
BestA = minimizer->GetParameter(0);
BestB = minimizer->GetParameter(1);
BestC = minimizer->GetParameter(2);

// Get the function value at the best fit.
chi2 = CHI2_lin_chiralfit_GP(BestA, BestB,BestC);

}

void lin_chiralfit_with1GP(double *Z, double **Z_CLUST, double *X, int NX, int Nclust,double &BestA,double &BestB,double &BestC,double &chi2)
{
	lin_chiralfit_withGP(Z,Z_CLUST,X,NX, Nclust, BestA, BestB, BestC, chi2);
}


void lin_chiralfit_with2GP(double *Z, double **Z_CLUST, double *X, int NX, int Nclust,double &BestA,double &BestB,double &BestC, double &BestD, double &chi2)
{

//Fit to ZETA=A+B*MPS^2+C/MPS^2
//MPS^2 propto (m1+m2)
// M=m1+m2 if POL
// M=MPS if MPS

// Reorganize data in global variables
NX_GLB=NX;
NCLUST_GLB=Nclust;
for (int ix=0; ix<NX_GLB; ix++)
{
        ZETA_GLB[ix]=Z[ix];
        X_GLB[ix]=X[ix];
	//printf("%8.5f %8.5f\n",X_GLB[ix],ZETA_GLB[ix]);
        for (int iclust=0; iclust<NCLUST_GLB+1; iclust++) ZETACLUST_GLB[ix][iclust]=Z_CLUST[ix][iclust];
}

// Minimize

TFitter* minimizer = new TFitter(2);
// MAKE IT QUIET!!

double p1 = -1;
minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);

// Tell the minimizer about the function to be minimzed
minimizer->SetFCN(minuitFunction_lin_chiralfit_2GP);

// Define the parameters
//arg1 – parameter number
//arg2 – parameter name
//arg3 – first guess at parameter value
//arg4 – estimated distance to minimum
//arg5, arg6 – ignore for now
minimizer->SetParameter(0,"A",0.0,0.01,0,0);
minimizer->SetParameter(1,"B",0.0,0.01,0.0,0);
minimizer->SetParameter(2,"C",0.0,0.01,0.0,0);
minimizer->SetParameter(3,"D",0.0,0.01,0.0,0);



// Run the simplex minimizer to get close to the minimum
minimizer->ExecuteCommand("SIMPLEX",0,0);

// Improvement fit doesn't work with double gp
// Run the migrad minimizer (an extended Powell's method) to improve the
// fit.
//minimizer->ExecuteCommand("MIGRAD",0,0);
//minimizer->ExecuteCommand("MINIMIZE",0,0);
//minimizer->ExecuteCommand("IMPROVE",0,0);
//minimizer->ExecuteCommand("MINIMIZE",0,0);
//minimizer->ExecuteCommand("MINOS",0,0);

// Get the best fit values
BestA = minimizer->GetParameter(0);
BestB = minimizer->GetParameter(1);
BestC = minimizer->GetParameter(2);
BestD = minimizer->GetParameter(3);

// Get the function value at the best fit.
chi2 = CHI2_lin_chiralfit_2GP(BestA, BestB,BestC, BestD);

}

double CHI2_lin_fit(double A, double B) {
        double CHI2=0;
        double F1;
        int NPAR=0;
        double *clust_aux=malloc(sizeof(double)*(NCLUST_GLB+1));
        for (int im=0; im<NX_GLB; im++)
        {
                NPAR=NPAR+1;
                for (int iclust=0; iclust<NCLUST_GLB+1; iclust++) clust_aux[iclust]=ZETACLUST_GLB[im][iclust];
                F1=(ZETA_GLB[im]-(A+B*X_GLB[im]))/error_jack_real(clust_aux,NCLUST_GLB);
                CHI2=CHI2+pow(F1,2);
        }

        if (NPAR<3)  {
		NPAR=3; //Avoid divide by 0
		printf("Warning: Divide by zero while computing chi2. Npar+1\n");
	}
        return CHI2/(NPAR-2);
//return CHI2;
}

void minuitFunction_lin_fit(int &nDim, double *gout, double &result, double *par, int flg) {
result = CHI2_lin_fit(par[0], par[1] );
}


void lin_fit(double *Z, double **Z_CLUST, double *X, int NX, int Nclust,double &BestA,double &BestB,double &chi2)
{

//Fit to ZETA=A+B*MPS^2+C/MPS^2
//MPS^2 propto (m1+m2)
// M=m1+m2 if POL
// M=MPS if MPS

// Reorganize data in global variables
NX_GLB=NX;
NCLUST_GLB=Nclust;
for (int ix=0; ix<NX_GLB; ix++)
{
        ZETA_GLB[ix]=Z[ix];
        X_GLB[ix]=X[ix];
        for (int iclust=0; iclust<NCLUST_GLB+1; iclust++) ZETACLUST_GLB[ix][iclust]=Z_CLUST[ix][iclust];
}
// Minimize
TFitter* minimizer = new TFitter(2);
// MAKE IT QUIET!!

double p1 = -1;
minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);

// Tell the minimizer about the function to be minimzed
minimizer->SetFCN(minuitFunction_lin_fit);

// Define the parameters
//arg1 – parameter number
//arg2 – parameter name
//arg3 – first guess at parameter value
//arg4 – estimated distance to minimum
//arg5, arg6 – ignore for now
minimizer->SetParameter(0,"A",1.5,0.001,0,0);
minimizer->SetParameter(1,"B",0.0,0.001,0.0,0);


// Run the simplex minimizer to get close to the minimum
minimizer->ExecuteCommand("SIMPLEX",0,0);
// Run the migrad minimizer (an extended Powell's method) to improve the
// fit.
minimizer->ExecuteCommand("MIGRAD",0,0);
// Get the best fit values
BestA = minimizer->GetParameter(0);
BestB = minimizer->GetParameter(1);

// Get the function value at the best fit.
chi2 = CHI2_lin_fit(BestA, BestB);

}


double CHI2_wprom_fit(double A) {
        double CHI2=0;
        double F1;
        int NPAR=0;
        double *clust_aux=malloc(sizeof(double)*(NCLUST_GLB+1));
        for (int im=0; im<NX_GLB; im++)
        {
                NPAR=NPAR+1;
                for (int iclust=0; iclust<NCLUST_GLB+1; iclust++) clust_aux[iclust]=ZETACLUST_GLB[im][iclust];
                F1=(ZETA_GLB[im]-(A))/error_jack_real(clust_aux,NCLUST_GLB);
                CHI2=CHI2+pow(F1,2);
        }

        if (NPAR<2) NPAR=2; //Avoid divide by 0
        return CHI2/(NPAR-1);
//return CHI2;
}

void minuitFunction_wprom_fit(int &nDim, double *gout, double &result, double *par, int flg) {
result = CHI2_wprom_fit(par[0] );
}

void wprom_fit(double *Z, double **Z_CLUST, double *X, int NX, int Nclust,double &BestA,double &chi2)
{

NX_GLB=NX;
NCLUST_GLB=Nclust;
for (int ix=0; ix<NX_GLB; ix++)
{
        ZETA_GLB[ix]=Z[ix];
        X_GLB[ix]=X[ix];
        for (int iclust=0; iclust<NCLUST_GLB+1; iclust++) ZETACLUST_GLB[ix][iclust]=Z_CLUST[ix][iclust];
}

// Minimize

TFitter* minimizer = new TFitter(2);
// MAKE IT QUIET!!

double p1 = -1;
minimizer->ExecuteCommand("SET PRINTOUT",&p1,1);

// Tell the minimizer about the function to be minimzed
minimizer->SetFCN(minuitFunction_lin_fit);

// Define the parameters
//arg1 – parameter number
//arg2 – parameter name
//arg3 – first guess at parameter value
//arg4 – estimated distance to minimum
//arg5, arg6 – ignore for now
minimizer->SetParameter(0,"A",1.5,0.001,0,0);


// Run the simplex minimizer to get close to the minimum
minimizer->ExecuteCommand("SIMPLEX",0,0);
// Run the migrad minimizer (an extended Powell's method) to improve the
// fit.
minimizer->ExecuteCommand("MIGRAD",0,0);
// Get the best fit values
BestA = minimizer->GetParameter(0);

// Get the function value at the best fit.
chi2 = CHI2_wprom_fit(BestA);

}

void wprom(double *Z, double **Z_CLUST, double *X, int NX, int Nclust,double &ZWPROM)
{

double *W=malloc(sizeof(double)*(NX));
double *aux_clust=malloc(sizeof(double)*(Nclust+1));
for (int ix=0; ix<NX; ix++)
{
	for (int iclust=0; iclust<Nclust+1; iclust++) aux_clust[iclust]=Z_CLUST[ix][iclust];
	W[ix]=pow(1/error_jack_real(aux_clust,Nclust),2);
}
ZWPROM=0;
double EWPROM=0;
for (int ix=0; ix<NX; ix++)
{
	ZWPROM=ZWPROM+Z[ix]*W[ix];
	EWPROM=EWPROM+W[ix];
}
ZWPROM=ZWPROM/EWPROM;

}

